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Abstract 

An ongoing problem in the study of a classical many-body system is the characterization of its 

^" 
O ■ equilibrium behaviour by theory or numerical simulation. For purely repulsive particles, locating 

O ■ 

the melting line in the pressure-temperature plane can be especially hard if the interparticle po- 
tential has a softened core or contains some adjustable parameters. A method is hereby presented 

-i— > ' 

c/2 ' 

that yields reliable melting-curve topologies with negligible computational effort. It is obtained 

by combining the Lindemann melting criterion with a description of the solid phase as an elastic 



continuum. A number of examples are given in order to illustrate the scope of the method and 
possible shortcomings. For a two-body repulsion of Gaussian shape, the outcome of the present 



approach compares favourably with the more accurate but also more computationally demanding 



£> , self-consistent harmonic approximation. 
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I. INTRODUCTION 

There exists as yet no comprehensive theoretical treatment of the solid-liquid phase tran- 
sition that can rival in quality and accuracy the smooth-cutoff formulation of the hierarchical 
reference theory of the liquid- vapour transition [2], which yields flat pressure vs. density 
isotherms in the coexistence region as well as distinct binodal and spinodal curves. Most 
theoretical and simulational strategies for detecting a point of solid-liquid coexistence invari- 
ably pass through the prior determination of the (Gibbs-) free energy for the two separate 
phases. Theoretical approaches include classical density- functional theories and thermody- 
namic perturbation theory J3j. In Monte Carlo simulations, the Frenkel-Ladd method and 
the Widom particle-insertion method may be employed, in conjunction with thermodynamic 
integration, in order to obtain accurate solid and liquid chemical potentials in the transition 
region |4|. On the far opposite side, lie a by now considerable number of empirical one-phase 
melting and freezing criteria which allow a rough estimate of the limit of stability for the 
given solid or liquid phase. Familiar examples are the Lindemann melting rule and the 
Hansen- Verlet freezing criterion, both relying on quantities that are computed numerically 
(the mean square displacement in the solid and the structure factor in the liquid). 

I hereby consider a semi-empirical method for the melting transition which, rather than 
being meant as a rule to provide a reliable estimate of the upper stability threshold of a solid 
with prescribed symmetry, is actually aimed at anticipating with little effort (i.e., without 
resorting to numerical simulation) at least the topology of the transition line, which may be 



useful especial 
are expected 



ly in those cases where multiple solid phases and/or reentrant-fluid anomalies 
5j. Assuming two-body forces between the particles, the idea is to treat the 
solid system as an elastic medium whose pressure-dependent moduli are determined at zero 
temperature from the potential. The upper limit of thermodynamic stability of the solid is 

n 

then taken in accordance with the Lindemann rule [6| . Aside from the approximate character 
of the Lindemann criterion, the main error in the estimate of the melting temperature T m (P) 
comes partly from neglecting all anharmonicities in the particle dynamics and partly from 
assuming the same elastic moduli at all temperatures. Both sources of approximation are 
expected to extend solid stability well beyond the actual threshold. In spite of this, the 
shape of T m (P) is reasonably well reproduced by this criterion, as I shall demonstrate for a 
number of model potentials (exceptions are anyway encountered - see below). 



After a reminder of elasticity theory in Section 2, I introduce the novel criterion of melting 
in Section 3 along with a few applications. In Section 4, I compare the indication of (what 
might be called) the elastic criterion of melting for the repulsive Gaussian potential with 
the outcomes of more refined approaches, based on the theory of the harmonic crystal and 
on the self-consistent harmonic approximation. Conclusions are postponed to Section 5. 

II. A BRIEF ACCOUNT OF LINEAR ELASTICITY 

The information gathered here is standard reference material which is preparatory to the 
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theoretical analysis that will be outlined in the next Section [8]. 

Consider the Bravais lattice {R„} of a crystal with N ^> 1 atoms (classical point particles) 
and let x ra = R n + u(R n ) be the actual position of the n-th atom (to simplify the notation, I 
choose Ri = in the following). In the simplest terms, the basic equations of linear elasticity 
stem from evaluating the total potential energy U of the crystal in the approximation where 
u(R) — u(R') is expanded to linear order in R — R'. For a homogeneous deformation, this 
amounts to replacing x na with 

3 

Rna + 2_^ UafiRnP (2. 1) 

0=1 

for a = 1, 2, 3, u a p being a constant. If the continuum limit is taken, we can think of u a p as 
just the constant value of du a jdx^. Concurrently, the deformation also modifies the crystal 
volume from Vo to V: 

— = det((5 Q/3 + u a p) = 1 + Y^ e *a + - ^2(e aa epp - e 2 ap ) + - ^ u 2 ap , (2.2) 

a a,/3 a,/3 

where third-order terms in the strains were neglected. In Eq. (12.2ft . 

e a /3 = ^{u a /3 + u f5a ) (2.3) 

is the strain tensor while u a /3 = (u a p — up a )/2. Our objective is to evaluate U up to 
second order in the u a p. Assuming a smooth spherically-symmetric pair potential <p(r) 
and specializing the analysis to crystals of cubic symmetry, a straightforward but tedious 
derivation yields: 

U 1 N 1 

u ~ ]v = 2 S ^(W) = u ° ~ p ( v ~ u °) + 2 v ° Xl (*** + e yy + ****) 

n=2 

+ ^oA 2 {e xx eyy + e xx e zz + e yy e zz ) + 2v X 3 (e 2 xy + e 2 xz + e 2 yz ) , (2.4) 



with v = V /N,u = (1/2) £n= 2 <K|Rn|), and 



1 N 

p = -^-Ei R «i^(i R ™i)- ( 2 - 5 ) 



6V ° „=2 



Moreover, three elastic constants (or Lame coefficients) appear in Eq. (12.41) : 

1 N X A 

Ai = - J P + ^ErRliO^lV / (|Rn|)-|R.|0'(|R.|)] ; 

i w \-2y2 

A 2 = ^+^ET^O R n|V / (|Rn|)-|Rn|0'(|Rn|)] ! 
•^n — * -tin 



J °n=2 



A 3 = A 2 -2P, (2.6) 

where X n , Y n , and Z n are the Cartesian components of R n . The three A's are the same 
quantities which are more commonly denoted Cn, Ci 2 , and C44, respectively. In Eq. (12. 4p . the 
term linear in the u a p and actually proportional to the trace of the strain tensor corresponds 
to the stress due to an applied pressure P. Equation ( 12. 51) links the lattice parameter (or 
the crystal volume Vo) with the pressure. The identification of P with the system pressure 
ensures consistency of Eq. ( 12.41) with the thermodynamic definition of pressure. 
A more general form of Eq. (12 .4p . valid for any temperature T, is the following: 

9 = 90 + 2^0 X] C «/37<5 e «/3 e 7<5 > ( 2 -7) 

where g is the Gibbs free energy per particle. Equation ( 12 .7p reduces to (12 .4p for T = and 
a crystal in the cubic system. The maximum number of independent elastic constants c Q/ g 7 ,5 
is 21 (taking Voigt symmetry into account), in fact they reduce to just three for crystals 
of cubic symmetry, five for crystals of hexagonal symmetry, and so on. For instance, for 
hexagonal solids Eq. (12. 7p takes the form 

g = g + 2v X 1 (e xx + e yy ) 2 + v \ 2 [{e xx - e yy f + 4e^] 

+ 2 V ^ e lz + 2v \ 4 (e xx + e yy )e zz + 4v A 5 [e\ z + e 2 yz ) , (2.8) 



with the following T = values of the Lame coefficients: 



1 N X* 

Al = T^ErRllO R nlV'(|Rn|)-|R„|0'(|Rn|)] 
i-^Vt) — * -rtn. 
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A -A P - 

A 2 " Ai--, 



1 7 4 

^ + 7T- E TrHI DR»I V(lRnl) - iRrMlRnl)] 
^0 o Ri 



n=2 

P 1 Aiiz 



^ ^2 7 2 



A 4 = v + ^ET^O R «lV'(|Rn|)-|Rn|0 / (|R„|)] ; 

A 5 = A 4 -P. (2.9) 

For tetragonal crystals, one similarly finds 

9 = 9o + y [Ai(eL + 4s/) + 2A 2e ra e ra + 4A 3 e^ 

+ A 4 eL + 2\ 5 (e xx + e yy )e zz + 4A 6 (e 2 xz + e 2 yz )] , (2.10) 

with zero-temperature Lame coefficients given by 

1 N X 4 

Al = -^H-^ErRliORnlVdRnD-lRnl^dRnl)] J 



n=2 



A 2 = P+^ETFf^ORnlV'(|Rn|)-|Rn|0'(|Rn|)] J 

^o n=2 |K n | 

A 3 = A 2 -2P; 

JL 7 4 
A 4 

'Zt'n 

n=2 

A 5 ' 



1 N 7 4 

"^ + y- £ TrHI R nl V(|Rn|) " |Rn|0'(|Rn|)] 
n=2 

? + ^tll [|Rn|V'(|Rn|) " |Rn|0'(|R 



^°„=2 



A 6 = A 5 -2P. (2.11) 

At T = 0, the expansion of the Helmholtz free energy P = Nf in powers of the strain- 
tensor components is the same as for U. For non-zero temperatures, the respective c a/ 3 7< 5 
are instead different (one thus distinguishes isothermal and adiabatic elastic constants). For 
any T, the Helmholtz free energy of a solid under arbitrary initial stress can otherwise be 
expanded to second order in the components of the displacement gradients u a p, 

= 2,S a pU a p + - \ SapjSUafiU-fS ', (2-12) 

^0 a * an 

a,p a,p,7,o 



alternatively, / can be written as a truncated power series of the Lagrangian strain param- 
eters, 

Vafi = r I U a p + Uf3 a + ^ U^gUrf J , ( 2 -13) 

with yet different coefficients in the linear and quadratic terms: 



VO 



a, 



a,/3,7,<5 

C a /3 and Sap^s 



(2-14) 
Ca^s + CpsSay Moreover, 



It is then a simple exercise to show that S a p 
for C a/ 3 = —P5 a/ 3, one finds that 

Ca/3-yS = Ca^S + P{5 a p^S — ^a-y^/38 ~ $a5$p~f) ■ (2-15) 

Equation (12. 15[) is useful for computing the A's at T > through numerical simulation since 
specific virial-like formulae exist for the C's [9|. 

An important issue is mechanical stability of a crystal phase, which is a prerequisite for 
its thermodynamic stability: an applied strain may destabilize the crystal, which in this 
case is really stable only at zero temperature. The elastic constants in Eq. (J2.7P must obey 
so-called stability conditions in order for the unstrained crystal to resist any infinitesimal 
deformation, i.e., in order for the crystal lattice {R„} to provide a minimum (not just an 
extremum) for g. Depending on the interparticle potential and on the pressure value, the 
crystal may or may not be mechanically stable, meaning that it does typically exist as a 
stable structure for T > only within one or more definite pressure ranges. 

Using Voigt symmetry, the elastic constants of a cubic crystal can be arranged in the 
6x6 matrix 



Ai A2 A2 


0^ 


A2 Ai A2 





A2 A2 Ai 








A 3 





A 3 


iy 


x 3 ) 



(2.16) 



,e 2 



and Eq. ( 12. 7ft becomes a quadratic form, g = g + (v o/2) J^ a b=1 c a be a eb with e\ = e a 
e yy e 3 = e zz,^4 = 2e yz ,es = 2e xz ,eQ = 2e xy . The eigenvalues (with multiplicities) of (12.1 6ft 
are A3 (3), Ai — A2 (2), and Ai + 2A2 (1), leading to three stability conditions: 

Ax + 2A 2 > ; A 3 > ; Ax - A 2 > . (2.17) 



The first two conditions amount to requiring the existence of the bulk and the shear modulus, 
respectively. The last inequality prescribes rigidity of the cubic solid against tetragonal 
shear. For a hexagonal crystal, a similar analysis yields four conditions, 

A 2 > ; A 5 > ; 8A1 + A 3 > ; X^ - \ 2 4 > , (2.18) 

becoming five for tetragonal crystals: 

A 3 > ; A 6 > ; Ai - A 2 > ; Ai + A 2 + A 4 > ; A 4 (Ai + A 2 ) - 2\ 2 5 > . (2.19) 

Tightly related to the subject of solid elasticity is the general harmonic theory of lattice 
dynamics. Consider a finite crystal with externally applied classical forces, and let the forces 
be restricted to the surface region so as to represent stresses applied to the crystal. Since 
the total force on each atom must vanish when the atoms are located at the equilibrium 
positions {R n }, the total energy at T = can be approximately written as 

[/ = (/() + ^EE $ ^ R - R'K(R)MR') (2-20) 

R,R' a,f3 

with Uq = C/(Ri, . . . , Rat), all anharmonicities being neglected. The $ coefficients in (I2.20p 
are second-order derivatives, 

^< R - R 'K aM^(R-) )o- (2 - 21) 

and, for a Bravais crystal, they are invariant under the exchange a •f-)- /3 because of the 
lattice inversion symmetry. Invariance of the energy value following a rigid translation of 
the crystal further leads to X^r^c^R) = for any a and (3. 
The equations of motion for the potential energy ( 12.201) read 

mu a (R) = -J2 <MR " R>/3( R , (2-22) 

R',/3 

where m is the particle mass, and are solved in terms of plane waves (the normal modes of 
vibration) , 

ea(q)e i[qR ^ (q)<] (a = 1,2, 3). (2.23) 

The N values of q lie within the first Brillouin zone (1BZ) of the lattice and are so chosen as 
to allow for the periodic repetition of the lattice outside its boundaries. Upon introducing 
the (real symmetric) dynamical matrix 

B a M) = J2 $ ^(R-)e 4q ' R = ~ Yl $ ^( R ) l l ~ «»(q ■ R )] > (2-24) 

R R 



the normal-mode amplitudes are found to obey the linear set of equations 

mu 2 (q)e a (q) = ^S a/3 (q)e^(q) . (2.25) 

For any q, the three eigenvalues of B a/3 (q), namely mu 2 (q) (s = 1,2,3), are real and we 
can always choose orthonormal eigenvectors, J2 a e sa(q)es'a(q) = <W- The explicit form of 
the dynamical- matrix components is B aa = r aa — T\ and B a p = T a/ 3 (a / /3), where 

i \ V^ ^'(|R«I) r-i /-dm 

T-i(q) = - 2^ ir> i t 1 - cos (<i ■ Rn)] ; 

^(^ = E TfTTT R «I V'(lRnl) " |R«|0'(|Rn|)] [1 - COs(q • R n )] . (2.26) 

A crystal dynamics is also associated with the approximation set by linear elasticity. It 
is drawn from the Lagrangian density (cf. Eq. (12.121) ) 

£ = 2^" 2 ( X ) ~~ zJ S <*pU a p - - 22 Safi-ySUapU-yS , {2.21) 

where p is the mass density. From Eq. (I2.27P one derives the equations of motion 

p^(x) = $> a/375 ^- , (2.28) 

whose solutions of are still plane waves with frequencies given by the secular equation 

det < ^ c a p 7 sqpqs - pu 2 (q)S ai > = . (2.29) 

In particular, one observes that the elastic waves are dispersionless, i.e., u 2 ex q 2 . For cubic 
crystals, the explicit form of Eq. (I2.29|) is: 

cnql + c ^(q y + ql) ~ pu 2 (C12 + c M )q x q y (c 12 + c M )q x q z 



(C12 + c u )q x q y c u q 2 + c Ai (q 2 x + q 2 ) - pu 2 (c 12 + c u )q y q z 

(c 12 + c u )q x q z (C12 + c u )q y q z c n q 2 + c^(q 2 x + q 2 y ) - pu" 2 



0. 



(2.30) 



III. A NEW MELTING CRITERION 

In this Section, a Gaussian field theory is formulated in order to describe the thermal 
properties of an elastic solid in the simplest possible terms. The aim is to obtain an approx- 
imate value for the mean square displacement (MSD) of crystal atoms that can be used to 
estimate the melting temperature of the crystal through the Lindemann criterion. 
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Consider for concreteness a crystal of cubic symmetry with iV = N1N2N3 atoms. Rather 
than assuming a homogeneous strain, I allow for a spatial dependence of atomic displace- 
ments and take the continuum limit. Then, the enthalpy H of the crystal at T = becomes 
(cf. Eg. (E1D): 



H = Hn 



1 

+ 2 
+ 2A 2 

+ A 3 



d 3 r { Ai 



v<> 



<9m 7 



9m r 9m 



dx 
9w r du 



+ 



9y 



+ 



du z 
dz 



y ^ ">x ^ "'z OUy OU z 



dx dy dx dz dy dz 



du T du. 



+ 



+ 



du r du. 



du„ du. 



dz 



dy 



(3.1) 



dy dx J V dz dx 
with the A's given by Eq. (12. 6p . Upon implementing periodic boundary conditions, the 
displacement vector is expanded in a series of plane waves: 

1 



u, 



£ 



M a (q)e iq ' r (conversely, w a (q) 



Vr 



d J nt a (r)e- ,qT 



(3.2) 



JV 



where, in terms of reciprocal-lattice primitive vectors, the wave vector q = Y2 a Icfaa with 
q a = m a /N a and m a = -N a /2 + 1, . . . , N a /2 (a = 1, 2, 3). Substitution of (H£2J into flSTTj) 
leads eventually to 

# = h + \v Q J2 E ^U(qR(q)^(q) (3.3) 



a,/3 



with 



4*/*(q) = [A 3 g 2 + (Ax - A 2 - 2X 3 )q 2 a ] 5 a$ + (A 2 + \ 3 )q a qp . (3.4) 

Next, I try to represent the thermal disordering of the crystal through a field theory where the 
basic variables are the M Q (r)'s and the statistical weight of field configurations is exp(—(3H). 
This choice is tantamount to the assumption of T-independent elastic constants, whose 
values are fixed at their (P-dependent) T = values. 

To compute the MSD, the following average is to be evaluated first: 
,~ ( ™ vv / m Vu* ga(qH(q) exp{-/3U £ q>0 E 7 ,^ ^(q)n 7 (q)^(q)} , , 

V°W«eW) /PMWexp{-/3UoE q >oE 7 ,^(q)^(q)^(q)} ' { ] 

where in both integrals the q's are restricted to half space (symbolically, q > 0) in order 
that {ReM a (q), IniM a (q)} can be treated as independent integration variables - namely, 
VuDu* = n q >o ri a d (ReM Q (q)) d (Im« a (q)). Using properties of complex- valued Gaussian 
integrals, one obtains 



<M Q (qK(q)) 



Vn 



(A- 1 ) 



a/3 



(3.6) 



Since the inverse of a symmetric matrix is also symmetric, the previous result actually applies 
for any q. Hence, the MSD reads 

1 / d 3 r M 2 (r)\ = ^^<K(q)| 2 > = ^^TrA- 1 (q). (3.7) 

V 0JV0 / 7 „ V ° q 

In the thermodynamic limit, the residual sum transforms into an integral over the 1BZ, 
which is more easily computed by replacing the zone with a (Debye) sphere of equal volume 
(the error committed is small), with the result: 



- \ dW(r)) = -^q D / d0 / d6 sin 6 J -f-f^ (3.8) 

VoJvo I k Jo Jo J2{9A) 



where qo = (67r 2 p) 1 / 3 and 



f x (0, 0) = A 3 (A 3 + 2Ai) + (Ai + A 2 )(Ai - A 2 - 2A 3 )(sin 4 0sin 2 cos 2 + sin 2 9 cos 2 9) ; 
/ 2 (0,0) = AiA^ + A3(Ai + A 2 )(Ai-A 2 -2A 3 )(sin 4 0sin 2 0cos 2 + sin 2 0cos 2 0) 

+ (Ai - A 2 - 2A 3 ) 2 (A! + 2A 2 + A 3 ) sin 4 9 cos 2 9 sin 2 cos 2 . (3.9) 

The parallel treatment for a harmonic crystal moves from 

U = Uo + f E E s ^(q)w«(q)^(q) , (3.10) 

q a,/3 

and eventually leads, through the same series of steps as before, to the following expression 
for the MSD, 

(^E M2 ( R )) = ^o£ D dgg 2 |d 2 fiTrfi- 1 (q), (3.11) 

which is more numerically demanding than (13. 8p because of the additional q integration 
present in (13.111) . 

We see from Eqs. (13. 8 p and (13.111) that the MSD increases linearly with T. According to 
the Lindemann criterion, the crystal melts when the MSD reaches a fraction L m w 0.1 of the 
nearest-neighbour distance a^^, from which the estimate of T m (P) follows directly. For face- 
centred cubic (fee), hexagonal close-packed (hep), and body-centred cubic (bec) crystals, the 
specific L m values are 0.15,0.10, and 0.18, respectively [lOj, lU|, while no systematic study 
of the typical values of the Lindemann ratio for other crystals has ever been undertaken, at 
least to my knowledge (hence, I assume L m = 0.1 indifferently for all such phases). If any of 
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the stability conditions is violated for a crystal under pressure Pq, then I take a zero melting 
temperature for the given solid at P = Pq. 

Universality of L m along the fluid-solid coexistence line is well established for fee, bec, 
and hep crystals. For other types of crystals no such information is available and this 
makes the T m estimated through what I shall call the elastic criterion of melting less reliable 
for these crystals. In general, the elastic constants get smaller and smaller on increasing 
temperature until they abruptly vanish on crossing the melting line. Hence, assuming the 
elastic constants to be independent of T is a major simplification that leads to systematically 
underestimating the MSD; moreover, also the neglect of anharmonic terms in the potential 
would likely contribute to enhancing the stability of the solid, with the effect that the T m (P) 
computed with the elastic criterion of melting will be larger than the actual value. One may 
reasonably expect that the extent to which the melting temperature is overestimated is 
roughly the same for all pressures so that at least the shape of T m (P) is got correctly. 

A first application of the elastic criterion is to the melting of the Lennard- Jones fluid, 
which is known to crystallize into a hep solid (unless the pressure is huge - larger than 800 
in reduced e/<r 3 units). For reduced pressures smaller than 20, the computed T m is a concave 
function of P, as expected [12j. For P = 1 and P = 10, the criterion predicts a melting 
temperature of 1.18 and 1.75, respectively, whereas the "exact" values from Ref. [12j are 0.78 
and 1.40. 

A more challenging test of the elastic criterion is offered by a recent simulation study |7J] 
of a system of particles repelling each other through the Yoshida-Kamakura (YK) potential, 



?yk(7J = eexp 



T N 


) 2 rn- 


(7/ 


° . 



o(l--) -6( 1 .-'-) \w'- (3.12) 



with a = 3.3. For reduced pressures smaller than 3, the phase diagram of the model is 
plotted in figure 3 of Ref. [7j . The same phase diagram but computed through the present 
melting criterion (with an enormous saving of time compared to simulation) is reported in 
Fig. 1. Here are shown the melting lines for a number of solid phases chosen among those 
stable at zero temperature. For each crystal, the melting curve is a single line or it consists 
of a number of disjoint pieces, one for each range of pressure/density where the stability 
conditions are met. It is worth stressing that the pressure range of mechanical stability 
of a phase is usually wider than the range of thermodynamic stability at T = 0, which is 
where the enthalpy of the phase is smaller than that of any other crystal phase. Hence, 
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FIG. 1: (Color online). Schematic phase diagram of the YK potential with a = 3.3 as drawn from 
the elastic criterion of melting. The melting lines of various solid phases are shown: fee (blue), bec 
(red), simple cubic (sc, black), sh (cyan), and /3-Sn (magenta). The dotted lines are the melting 
curves for the fee and bec crystals as derived from the harmonic approximation, see Eq. (|3.1ip . In 
the inset (top panel), a comparison is made with the exact coexistence boundaries of the model 
(black dots and thick solid lines) [7|. From low to high pressure, the stable phases up to P = 3 are 
fee, bec, and /3-Sn. 

the stability boundaries dictated by the elastic criterion do not generally coincide with the 
actual thermodynamic thresholds. 

On approaching a stability boundary, the MSD of Eq. ( 13. 8 p blows up and the melting 
temperature drops continuously to zero. The line of fluid-solid coexistence would correspond 
to the upper envelope of the melting curves for the various solids. It is clear from Fig. 1 that 
the gross features of the phase diagram of the YK fluid are well reproduced by the elastic 
criterion, the main error being in the regular overestimation of the melting temperature. 
The greater stability of the /3-Sn phase over the simple hexagonal (sh) solid in the pressure 
range between roughly 3 and 7 might be just accidental, related to the choice of the same L m 
for both. The harmonic approximation works quantitatively better (since at variance with 
linear elasticity no large-wavelength limit is implied) but it takes a much longer computer 
time to calculate the MSD. 
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0.05 



FIG. 2: Low-pressure phase diagram of the YK potential with a = 3.3 according to the theory de- 
tailed in the Appendix. Solid-solid coexistence points are depicted as small dots, whereas triangles 
and the full square are solid-fluid coexistence points. When there are more than one crystal phase 
of a given type, a Roman numeral distinguishes between them (e.g. /3-Sn-I and /3-Sn-II; the second 
fee phase is stable for pressures out of the range shown). The dotted lines through the open dots 
are the coexistence loci of the model from Ref. 171. 



To better appreciate the quality of the elastic criterion of melting, it is worth considering 
what would be the phase diagram of the YK potential with a = 3.3 according to a theory of 
fluid-solid coexistence based on the use of the cell-theory approximation for the solid and the 
Mansoori-Canfield theory for the fluid (see the details in the Appendix). We see from Fig. 2 
that this theory predicts a direct transition from bec to sh at high temperature, a possibility 
which was not actually considered in the simulation; however, the melting temperature of 
the YK fluid is overestimated by the theory to roughly the same extent (~ 100%) as it is by 
the elastic criterion, a fact that alone casts some shadows on the reliability of the theoretical 
phase diagram. 

It is instructive to look at the shape of some representative phonon branches of the bec 
crystal of YK particles for p = 0.6607 (P ~ 2.76), i.e., where the bec solid is about to 
become unstable at zero temperature owing to the fact that C44 is almost zero and actually 
negative for larger pressures. This instability is caused by phonon softening at the T point: 
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FIG. 3: Yoshida-Kamakura potential (|3.12p with a = 3.3: phonon branches of the BCC crystal 
for p = 0.6607 along a number of high-symmetry lines in q space. Along the TN path, one of the 
branches is seen to soften at the T point due to the vanishing of c 44 . 

along the path from T to N, one of the phonon branches satisfies mu; 2 (q) ~ (c 44 /p)(g 2 + q^) 
for q — y (see Fig. 3). 

Upon varying the value of a in Eq. f 13 .12)) . one can follow the evolution of the YK phase 
diagram through the elastic criterion of melting [7| . For large a values, the inverse-power- 
fluid limit is recovered; for a ~ 7, there appears a region of bcc stability between the low- 
and high-density fee solids; on decreasing a more and more, the stable-bcc region gradually 
shrinks until, for a ~ 4, a gap opens between the bcc and high-density fee regions, signalling 
the stabilization for intermediate pressures of one or more crystals of symmetry other than 
cubic. The opening of the gap is preceded by the onset of reentrant melting, which first 
occurs for a « 5. 

Another instance of core-softened repulsion is provided by the modified inverse-power 
(MIP) potential studied in Ref. 131 ] . The following one-parameter family of potentials is 
being considered: 



/cr\ n (' r ) 
(f*Mip(r) = e ^-J with n(r) 



12 <^ 1 



aexp 



-5(1-- 
a 



(3.13) 



where < a < 1 is a softness parameter, i.e., a number fixing the extent to which the 
inverse-power exponent deviates from 12 in the close neighbourhood of a. Upon increasing 
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FIG. 4: (Color online). Modified inverse-power potential with a = 0.8. Left: Numerically- 
computed phase diagram (reprinted from Ref. [13]; the dots are melting points as obtained by 
the heat-until-it-melts method while the vertical dotted lines are putative solid-solid boundaries 
as extrapolated from exact total-energy calculations at T = 0); right: same phase diagram as 
predicted through the elastic criterion of melting (blue, fee; red, bec; black, sc; cyan, sh; magenta, 
/3-Sn; the blue and red dotted lines are the melting curves for the fee and bee crystals, respectively, 
as drawn from the harmonic approximation). 

a, the potential core softens more and more, with the effect of destabilizing both the fee 
and the bec order for intermediate densities. This is accompanied by reentrant melting and 
by the appearance of one or more low-coordinated crystal phases in the pressure gap left 
open by bec and fee. In the left panel of Fig. 4, I report the phase diagram of the MIP 
fluid for a = 0.8 as obtained from Monte Carlo simulation through the heat-until-it-melts 
method [13j; the same melting lines but derived from the elastic criterion are plotted in 
the right panel of Fig. 4. Again, we see more than one correspondence between the present 
melting criterion and the simulation results. 

However, there are also instances (arguably not so common) where the elastic criterion 
fails badly. This occurs when a crystal that is predicted by linear elasticity to be unstable 
at T = is in fact stabilized in a range of temperatures, somewhat counterintuitively, by 
virtue of anharmonic effects. In a case of these, anharmonicity manages to make a crystal 
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TABLE I: MIP potential for a = 0.6, elastic constants of the BCC crystal at the reduced densities 
p = 0.7 and p = 1. The exact T = values derived from Eqs. (|2.6[) are compared with their MC 
estimates at T = 0.001 (for samples of iV = 686 particles and equilibrium trajectories of as many 
as 2 x 10 5 MC moves per particle). While the BCC crystal would be mechanically unstable at 
p = 1 according to elasticity theory, it is actually found perfectly rigid to thermal fluctuations in 
numerical simulation owing to the stabilizing effect of the anharmonicities in the potential. 



p = 0.7 


p = \ 




Cu C12 C44 


Cll C12 


c 44 


T = 10.06156 8.64219 2.94964 
MC 10.046(1) 8.632(1) 2.937(1) 


18.92754 13.80805 
18.1(3) 13.89(3) 


-0.41631 
0.955(4) 



phase rigid to small deformations in spite of the violation of the stability conditions of 
elasticity. I found one case of these for the MIP potential. The T = calculation of the bcc 
elastic moduli for a = 0.6 predicts a gap of stability in the density range from p = 0.910 
(P ~ 5.761) to p = 1.066 (P ~ 8.168), whereas for e.g. p = 1 (P ~ 7.021 at T = 0} 



el). What is happening then? 
14j : strong enough anharmonic 



Monte Carlo simulation clearly indicates that the bcc solid is stable up to T ~ 0.105 [13] 
(all quantities in reduced units). A numerical calculation of the elastic constants for p = 1 at 
very low temperature (T = 0.001) with the method of Ref. [9J] indeed reveals large deviations 
from the T = values, which is not the case for e.g. p = 0.7 (P ~ 2.846 at T = 0), where 
the agreement with linear elasticity is much better (see Tab 
The similar situation with Calcium sc phase provides a clue 
terms in the Hamiltonian (classical or quantum) may succeed to convert imaginary phonon 
frequencies into real ones, thus allowing the alleged unstable solid to become mechanically 
(and thermodynamically) stable. 

IV. THE MELTING CURVE OF THE GAUSSIAN-CORE MODEL 

The Gaussian-core model (GCM) fluid (i.e., classical point particles interacting through 
a repulsive Gaussian potential in three dimensions) gives the opportunity to compare the 
relative efficacy of various empirical melting rules, all rooted in the use of the Lindemann 
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criterion. In particular, we shall figure out the merits and drawbacks of the self-consistent 



harmonic approximation (SCHA) [15|, which for many years represented a popular theoret- 



ical alternative to exact free-energy calculations. 



Besides a fluid phase, the GCM shows two distinct, fee and bec solid phases 16J. At 
T = 0, the fee solid transforms to bec for P = 0.05529. At higher pressures and for T > 0, 
the bec solid undergoes reentrant melting: T m (P) is an increasing function for P < 0.136 
while being decreasing otherwise, further vanishing in the limit of infinite pressure. The fee 
and bec melting lines as predicted by the elastic criterion are reported in Fig. 5, together 
with those obtained from the harmonic approximation. In the same picture, the outcome 



of a variational treatment [17| and the numerically-computed coexistence lines [18( are also 
plotted for comparison. Clearly, the simple elastic criterion is able to account for the main 
characteristics of GCM melting, though the fee and bec melting temperatures are again 
found to be about twice larger than the actual values and the threshold where the fee solid 
is overcome in stability by the bec phase remains vague, much overestimated by the puta- 
tive fee reentrant-melting line. Quantitatively speaking, the harmonic approximation and, 
especially, the variational theory provide more valid alternatives to free-energy calculations. 

The SCHA is a theory for the thermal attenuation of phonon energies that aims at 
introducing elements of anharmonicity in an otherwise harmonic set-up. It provides an 
internal, self-consistent condition for its own validity which had sometimes been interpreted 
as an indication of the maximum temperature at which the crystal can be superheated. 
When used in combination with the Lindemann rule, the SCHA provides an independent 
melting criterion. Before illustrating the specific prediction for the GCM, I present a brief 
introduction to the SCHA. 

The formal justification of the SCHA lies in the use of the variational method of statistical 
mechanics. The strategy is focussed on determining the "optimal" harmonic approximation 
to the real Hamiltonian at the given temperature T, which is generally not its harmonic 
part. The crucial assumption is that of an integrable pair potential </>(r), endowed with 
a Fourier transform 0(q). This automatically excludes hard-core potentials, for which the 
SCHA theory cannot be formulated. The average of the system potential energy over a 
reference harmonic system U^m, having the same potential-energy minimum as the system 
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FIG. 5: (Color online). The Gaussian-core model phase diagram as determined through vari- 
ous methods: exact free-energy calculations (solid black lines); variational method (dotted black 
lines); elastic criterion of melting (long-dashed blue and red lines - blue, fee; red, bec); harmonic 
approximation (dashed blue and red lines). 



of interest but different phonon frequencies u s (k) and normal-mode amplitudes e s (k), is 



(U) 



Vq 



harm 



^E 



i- -j 



d3g T 'q)e iq -( R *- R ^V^ (q ' (Ul ~ Uj))2 )i 



(27T) 



(4.1) 



where the prime over the sum means i ^ j and 

harm AT 



1 <( q ■ (u, - Uj )) 2 > harm = k -^- £ (q ■ e *( k )) 2 - ^m 3) s D(q ' {R}) • (4 - 2) 



The best approximation to the Helmholtz free energy of the system within all conceivable 
harmonic interactions is given by the minimum of the Gibbs-Bogoliubov functional, 



3 
2' 



F[Hhaxm] = i*harm + (H — -ffharm)harm ~~ -^harm + (^)harm — ^0 ^ T^fcgT , 



where the Helmholtz free energy of the reference system reads 



(4.3) 



^harm = U + 3Nk B T\n \ -^ ] + — - ^ In 



,V3 



rcksT 



(4.4) 
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with A the thermal wavelength. Using the frequencies w a (k) as variational parameters, they 
are eventually obtained as the solutions of the SCHA equations 

™ s 2 (k) =v J2 (e- 1 ^ ~ 1) / ^03 (q ■ e s (k)) 2 0(q) e -*i^ e -^W) . (4.5) 

In practice, the target temperature T is reached in steps, where at every step of the calcu- 
lation the equations (J4.5P are solved iteratively until the left-hand side equates to a certain 
degree of precision the right-hand side. At low temperature, a good starting point of the 
iteration are the system own frequencies. Observe that, thanks to symmetry considerations, 
a (congruous) number of k vectors in a small fraction of the 1BZ will suffice for the calcula- 



tion of a sum like that in D (e.g. just 1/48 of the full 1BZ for the FCC lattice) 19]. Once 
D is obtained, the matrix 

/d 3 o ~ 

77^ <?«9/30(q) cos(q • Rj)e~ D (4.6) 

is diagonalized in order to extract its eigenvalues mw s 2 (k) and eigenvectors e s (k), and this 
completes a single iteration step. 

The main limitation in the use of the SCHA method is computational, due to the necessity 
of solving numerically a large number of times the integral in ( 14. 6 p to a high degree of 
precision. In the GCM case this integral can be computed analytically and this enormously 
speeds up the whole procedure. Even in this favourable situation, computing a single melting 
point by the SCHA method takes a time typically three orders of magnitude larger than if we 
apply the elastic criterion, which performs the calculation in a few hundredths of a second on 
a fast PC. In general, for a given density p the self-consistent calculation of the frequencies 
w s (k) and the respective MSD can be accomplished only up to a certain temperature T,(p), 
which is called the instability temperature. Beyond this temperature, no self-consistent 
solution of the Eqs. (I4.5P is found. Moreover, depending on p the ratio of the MSD at Tj 
to oat at may even exceed L m (in this event, I assume T m = %). The SCHA results for the 
GCM are reported in Fig. 6. Compared to the outcome of the elastic criterion, the SCHA 
estimate of the GCM melting temperature is better for all low to intermediate densities; the 
SCHA is instead unable to reproduce the large-density tail of T m (P) since beyond a density 
of 0.57 I find no self-consistent solution of the Eqs. (14. 5 p and the relative T m hence drops to 
zero. It is worth adding a final remark about the SCHA instability thresholds for the GCM 
model. As we see from Fig. 6, the Tj for fee is higher than it is for bec, in sharp contrast with 
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FIG. 6: (Color online). Gaussian-core model phase diagram: the outcome of the SCHA (blue and 
red dots, joined by solid straight lines - blue, fee; red, bec) is compared with that of the elastic 
criterion of melting (long-dashed blue and red lines). The SCHA instability temperatures for the 
fee and bec solids are also plotted as blue and red crosses, respectively, joined by dotted straight 
lines. Finally, the solid black lines mark the exact coexistence loci. 

the sequence of melting thresholds. In fact, the SCHA instability at, say, pa 3 = 0.2 occurs, 
for both phases, where the root mean square displacement (rmsd) for the reference system 
is roughly a fraction 0.23 of the nominal NN distance (rjvjv)- But the rate of growth with 
temperature of rmsd/rjvTv is slightly larger for bec, with a pronounced acceleration above a 
level of about 0.17 for bec and 0.20 for fee; hence, the rmsd/ 'r^N of the bec crystal reaches 
the values 0.18 (melting) and 0.23 (instability) both within the range comprised between 
the fee melting and instability temperatures. 

V. CONCLUSIONS 

Through the use of representative model potentials, I managed to show that a simple 
melting criterion based on the Lindemann rule and a description of the solid as an elastic 
medium is able to capture, with negligible computational effort, the overall characteristics 
of the system melting line. In more quantitative terms, the criterion overestimates the 
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melting temperature by roughly a factor of two for fee and bee solids, independently of the 
pressure value. For other crystals, the prediction of the criterion is less reliable, mainly 
due to the uncertainty on the value of the Lindemann parameter and its actual pressure 
dependence. In fact, the value of the elastic criterion of melting is more of a heuristic kind, 
i.e., of guidance for fastly detecting the existence of anomalies in the melting line, as in 
case of reentrant-fluid behaviour in a system where the softness of the particle core can be 
made to vary by tuning an appropriate parameter in the potential. The accuracy obtained 
by the elastic criterion in predicting the overall appearance of the phase diagram can be 
comparable to that of more sophisticated (two-phase) theories of fluid-solid coexistence, 
as I showed for one instance of core-softened interaction. For a Gaussian repulsive core, 
I compared the outcome of the elastic criterion with the harmonic approximation, as well 
as with the more effective but also more numerically demanding self-consistent harmonic 
approximation. Though moving upward in the hierarchy of theories generally improves the 
estimate of the melting temperature for all pressures, the gain in accuracy is only marginal 
and, more important, the topology of the melting line stays unaltered. Hence, at least for 
the Gaussian potential, a description in terms of zero-temperature elastic constants is by far 
sufficient to anticipate the essential features of the melting behaviour and no better theory 
is strictly necessary. 

Appendix A: A statistical theory of the fluid-solid transition 

In this Appendix, a theory of fluid-solid coexistence is formulated for purely repulsive 
potentials, where the fluid phase is described through the variational approach by Mansoori 



and Canfield 20| while the statistical properties of each solid phase are modelled through a 
cell theory. 

Assuming the hard-sphere (HS) fluid as reference, it derives from the Gibbs-Bogoliubov 
inequality that the exact Helmholtz free energy per particle / of a system with potential 
<f>(r) is bounded from above by 

/* (T, v- a HS ) = /hs + | / d 3 r g H s(r)w(r) with w(r) = cj>{r) - HS (r) , (A.1) 

where p = 1/v is the number density and gns{ r ) is the HS radial distribution function 
(RDF). In Eq. flA.lj) . the HS particle diameter <ths is left unspecified; the best approximant 
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to / is obtained by minimizing (1A.1J) with respect to cths- Although the HS equation of state 
is not known exactly, a good approximation is the Carnahan-Starling form [3] from which 
the HS free energy follows as 



/ HS = k B T [ln(pA 3 ) - 1] + k B T ^_^ 



(A.2) 



with 1] = (7t/6)p(Ths- To obtain an estimate of the HS RDF, one resorts to the Percus-Yevick 
approximation |3J since then the direct correlation function Cns( r ) — Co(r/o"Hs; v) i s known 



in a closed form: 

Cq(x) 



— Ao — Aix — A3X 3 , x < 1 
, x> 1 

(1 + 2 V ) 2 



with A, 



Ai 






(I-77) 4 ' - 1 "' (1-77) 

The Ornstein-Zernike relation then yields gns{ r ) — 5 , o( r / cr Hs! v) with 



flbW 



1 + 



7T 



, sin(A;x) co(fc) 



d k k 2 

n /ex 1 — 24r/co(A;) 



and Co (A;) 



dxx' 



, sin(fcc} 
kx 



(A.3) 



c (x) . (A.4) 



The variational free energy (IA.1I) can then be written as 



/* = k B T [ln(pA 3 ) - 1] + fc^T ^- *$ 4 j 2// / d,-,- 2 , /n (,-: ,,)o(.ro us ) 



7^ 



A 



3A; B rin - + k B T [ln(pcx 3 ) - l] + Af* 



a 



(A.5) 



where a is an arbitrary length unit. Called <7ns{T, v) the optimal <7hs value and observing 
that Af* depends on v only through r\ (i.e., Af*(T, v\ cths) = <f(T, 77(11, cr H s); chs)) 5 the fluid 
chemical potential can be approximated as p = f + Pv, where / = f*(T,v; aus(T,v)) and 
P = —df/dv. In order to calculate P, one considers that 



dAp 



da 



As a result, 



HS 



Pv = — f 



T,ti 



dv 



whence -7— 
ar\ 



CHS CV 



r.^HS 



T,», 



knT — V 



dAf 



dv 



k B T- 



3r? da HS 
er H s <9<p 



(A.6) 



T 



3 cV H s 



T,»J 



A; B T - — -^ 

3 f ./! 



dxx 3 5'o(a;;r7)0 / (a;crHs) 



(A.7) 



This completes the derivation of an approximate expression of the fluid chemical potential 
to be compared with the chemical potential of the solid phase. 
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FIG. 7: Left: phase diagram of the YK potential with a = 2.1 according to the theory detailed 
in the Appendix (for the meaning of Roman numerals, see Fig. 2 caption). Right: zoom on the 

dots, whereas triangles, 



23] for a comparison with 



low-pressure region. Solid-solid coexistence points are depicted as sma 

squares, tripods, and crosses are solid-fluid coexistence points. See Ref. 

the prediction from Monte Carlo simulation. The open dots with error bars give the location of 

number-density maxima within the Mansoori-Canfield description of the fluid phase. 

As far as the solid sector of the phase diagram is considered, I first determine the stable 
phases at T = through a series of total-energy calculations for a large number of candidate 
crystal structures (see Ref. [2l| for more details). To obtain a rough estimate of the crystal 
chemical potential at T > 0, I use the simple Lennard- Jones- Devonshire cell theory [22J. In 
this theory, a crystal partition function of effectively independent particles is written down 
where any given particle, which can be found anywhere in its own Wigner-Seitz cell (WSC), 
is acted upon by the force exerted by the other N — 1 particles, placed at equilibrium lattice 
positions. In practice, the canonical partition function of a crystal is approximated as 



where 



Z = W / d3ri • ■ ■ / ^ ex P I - E ^d/{k B T) 1 

A JWSCi ^WSCiv [ j J 

(r) = \ J2 0(|Ri - Ril) + E ^l Rl + r - R iD - <KI R i " R;l 



(A.8) 



(A.9) 



m 
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Taking 
$(r) = ^0(|R 1 + r-R j |) and *(r) = -~ J] l R i + r " *W(|R<i + r - R-il) , (A. 10) 



3 



a direct calculation offers 



F 
N 



and 



-fc B Tln^ + -$(0) withv/= / d 3 r exp{-/3[$(r) -$(0)]} (A.ll) 

A 2 iwsc 

fI = L + Pv = 3k B T\n - - k B T fin ^ - l) + ~ [$(0) + *(0)] 

J wsc d 3 r [g(r) - vfr(O)] exp {-/3[$(r) - $(0)]} 

j wsc d3rexp{-/?[$(r)-<l>(0)]} ' l ' j 

Fig. 7 (left panel) shows the phase diagram of the Yoshida-Kamakura potential (13.121) for 
a = 2.1 as mapped out in the way just explained. A zoom on the low-pressure region of 
the phase diagram is presented in the right panel of Fig. 7. Compared to the exact phase 
diagram of Ref . 23|] , we see that the theory correctly accounts for the succession and extent 
of solid phases (with the unique omission of the cI16 solid), though still overestimating the 
values of the melting temperature by approximately 100% for all pressures. In the same 
picture, I also plotted the line encompassing the region of density anomaly as computed 
within the Mansoori-Canfield theory. The shape of this line compares well with that of the 
same line as obtained from simulation. 
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